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Reformulating hyperdynamics without using a transition state theory (TST) dividing surface makes it 
possible to accelerate conventional molecular dynamics (MD) simulation using a broader range of bias 
potentials. A new scheme to calculate the boost factor is also introduced that makes the hyperdynamics 
method more accurate and efficient. Novel bias potentials using the hyper-distance and the potential 
energy slope and curvature along the direction vector from a minimum to a current position can 
significantly reduce the computational overhead required. Results simulating an atomic force microscope 
(AFM) system validate the new methodology. 

I. INTRODUCTION 

When a dynamical system is confined in a potential energy basin separated by large energy 
barriers (» ^^7), the system stays in this basin for a very long time compared to a typical atomic 
vibrational period before hopping to other basins. If the potential energy has multiple such minima 
(metastable states), the system evolves through inirequent transitions irom one metastable state to another. 
In physical systems dominated by inirequent events, information about the waiting times at each state, the 
transition mechanisms leading to other states, and their relative probabilities is essential to understand and 
predict their behavior. 

Direct dynamics simulation methods like molecular dynamics (MD) ' do not assume any prior 
knowledge about how a system will evolve in time and can therefore be used to simulate as yet unknown 
transition mechanisms. However, since the overall time scale accessible by the conventional MD 
simulation method is restricted by the atomic vibrational period due to issues of numerical stability, it has 
been difficult to simulate these inirequent event problems using the MD methodology. Alternatively, if 
we were able to enumerate every transition mechanism and its rate for all the states the system visits 
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during its evolution, this would allow us to advance the system from one state to another without 
following detailed trajectories in configuration space. This is the fiindamental assumption of Kinetic 
Monte Carlo (KMC).^ Moreover, when two adjacent minima are known, the transition rate between these 
two states can be calculated using transition state theory (TST) ^"^ if a proper dividing surface that the 
system crosses in transitions can be constructed. However, determining all transition mechanisms 
becomes more and more intractable as system complexity increases. Therefore, performing dynamics 
simulations on time scales reachable by KMC has long been a goal. 

In recent years several novel methods to extend the MD time scale have been proposed. These 
methods include hyperdynamics,^' the parallel replica method," and temperature-accelerated dynamics 
(TAD).'^ In the parallel replica method multiple replicas of a given system are simultaneously simulated 
on different processors so that the transitions can be accelerated up to a factor corresponding to the 
number of processors." TAD uses a higher temperature MD simulation to more effectively detect 
transition pathways at a given lower temperature, but it assumes that the harmonic TST holds and 
requires the saddle point to be found for each transition pathway. A more detailed review of these 
methods is also found in Ref. 14. 

Taking into account both the achievable boost factor and the degree of approximation, 
hyperdynamics is perhaps the most attractive acceleration method. In hyperdynamics, ' a given potential 
energy fiinction is modified such that the energy barriers are reduced while the characteristic dynamics 
are preserved. In principle hyperdynamics simulation can advance the system at an accelerated pace while 
preserving the correct relative transition probabilities under the assumption that the TST transition rates 
are equivalent to the actual rates. Furthermore, the acceleration rate can be calculated concurrently during 
the simulation. However, constructing a computationally efficient bias potential for a hyperdynamics 
simulation can be difficult. In Voter's original bias potentials the lowest eigenvalue and the corresponding 
eigenvector of the Hessian matrix are used, but calculating the eigenvalue and its derivative require 
significant computational overhead. To avoid these excessive computations several simplified bias 
potentials have been proposed, and recently Miron and Fichthom proposed a bias potential called the 
bond-boost method using the bond length changes to detect a transition without significant computational 
overhead.'^ The bond-boost method utilizes the characteristic of bond-breaking that most solid-state 
systems undergo when making transitions to construct a bias potential. However, this method gives rise to 
a significant force discontinuity due to the envelope function introduced to enforce that the bias potential 
is zero at the dividing surface. It can also introduce spurious energetic minima within an individual bond. 

Although several bias potentials have been proposed thus far, all either have significant 
computational overhead, which degrades the achieved boost factor, or lack the generality required to 
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detect a diverse range of transitions. Thus, finding a bias potential that provides a large boost without 
significant overhead remains a challenging problem. 

In this manuscript we begin by reviewing the fiindamentals of the hyperdynamics method and 
then proceed to reformulate the method in a rigorous way that obviates the need to construct a TST 
dividing surface. Special emphasis is put on the importance of accurately calculating the boost factor. 
This will set the stage for devising new bias potentials using local variables in place of or in addition to 
the lowest eigenvalue of the Hessian matrix. Finally, the new methodology is tested with an atomic force 
micro scope (AFM) system modeled using the Lennard- Jones potential. 



II. REVIEW 



A. Rare events and transition rate 



Let us consider a system in an ensemble with constant boundary conditions such as the canonical 
ensemble (NVT). The characteristics of this system can be described by the Hamiltonian i/ defined by 

H(r,v)^V(r) + Kiv) , (1) 
where r is the 3iV-dimensional position vector in the configuration space (N is the total number of 
particles), v is the 3A^-dimensional velocity vector. Vis the potential energy, and A" is the kinetic energy. 
The probability density distribution p{r, v) in the phase space is given by 

p(r, v) = ^Qxp(-^H(r, v)) , (2) 

where Z is the partition function defined hy jdfjdv e~^" and j3 = l/kgT; ks is the Boltzmann 
constant and T is the temperature. Then, the ensemble average of any observable 0(r, v)is expressed as 

{0{r,v)) = ^d7^dv 0{r,v)p{r,v) = ^^d7^dv 0{r,v)e~'^" . (3) 

The dynamics of this system {r{t),v{t)} , a trajectory in the phase space satisfying the probability 
density in Eq. (2), can be modeled, for example, by coupling an isolated system to the Nose-Hoover 
thermostat.'*"'^ 

If the potential energy has multiple local minima separated by high energy barriers (» ^^7), each 
local minimum in configuration space and its neighborhood, around which the probabilities are 
concentrated, defines a metastable state and the characteristic dynamics of the system consists of 
infi-equent transitions from one metastable state to another as illustrated in Fig. 1. If the waiting time at 
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each state is comparable to the observation time, the time average of an observable is not given by Eq. (3), 
but depends on the specific transition paths the system follows. We assume that the metastable states each 
have a corresponding non-overlapping configuration space volume, called a state voliune, so that they are 
distinguishable from each other and do not share any position vector as illustrated in Fig. 1. 

Now we shall describe how to calculate the waiting times and the relative transition probabilities. 
Let us consider a transition from a state A to one of two neighboring states B and C (hereafter we use 
capital letters to refer to either a state or a state volume unless an ambiguity would arise). Then, a 
transition from state A occurs only when the system trajectory enters the state volume B or C after leaving 
the state volume A rather than re-entering itself The system can cross and recross the boundary of the 
state volume A several times before making a transition (see Fig. 1). Thus, when we say that the system 
enters, leaves, and stays at state A, we refer to the initial entrance event to the state volume A, the final 
leaving event, and the time period between these two events, respectively. This definition of a transition is 
essentially the same as that in Ref. 8. 

The waiting time t^is defined as the total period of time since the trajectory initially enters the 
state volume A until the trajectory finally leaves the state voliune A. If the system leaves the state volume 
A only when making a transition, the waiting time is the same as the time period the system spends inside 
of A. However, because in practice it is difficult to construct such an ideal state volume, we consider the 
possible recrossings of the boundary of A before the trajectory finally leaves it. Then, can be expressed 
as 



where t^^ is the total time the system spends inside of A and „ is the total time the system spends 
outside of A before making a transition. The totality of the positions outside of A where the system can 
visit before making a transition is referred to as the transition region of A and is denoted as A . Note that 
if the state volume A is well chosen t^j»t^g . 

Since in infrequent events the waiting times are long and uncorrelated we can assume that t ^ has 
a Poisson distribution expressed as ^ 



where P (• • •) is a probability density function, and Ra is a rate constant that characterizes the transition. 
We can also show that the following relation holds. 



(4) 



(5) 




(6) 
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where E (• • •) is an expectation value. If we can construct an infinite trajectory such that the system visits 
and escapes fi"om state A a large number of times, the average waiting time can be calculated by 



r^=lim^ = E(?J , (7) 

where t is the total time elapsed by the system; N^^ is the total number of the transitions from state A in 

this time interval; is the waiting time during the Ath stay at state A. Hereafter we use bars to refer to 
the average of samples obtained in this infinite trajectory as in Eq. (7). By rearranging Eq. (7), we obtain 



_ ^ _ Pa 



?^=lim ^=^ , (8) 

where 



r— >oo 7" AT II 



u^^^lim^ , (9) 



k=\ 



= lim ^ — . (10) 

where v^^ is the mean frequency of the transition from state A and is the ratio of the total waiting 
times at state A to the entire time interval, can be expressed as 

Na^ Na^ Na^ 

p^^]im^^^]im^ '-^ ^p^j+p^^, , (11) 

where 

Na^ N, 



k 

A,o 



p^^^ and ^lirn^^^^ . (12) 

By definition, the state volume A does not overlap with other state volumes so that in an ergodic system 
we have the following relation, 



r e ^ 



However, in general A can overlap with B or C so that we have 
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p^^,=^^ < or p^^^=cy-^—— (c^<l). (14) 



r— >co 77 I rl V ^ ■ ^ ^ T Q 



If A does not overlap with other transition regions or the probabilities of the overlapping regions are 
negligible, then « 1 . In the derivations belovv' v^^e ignore the possibility that c^^^X. The procedures to 
calculate the mean frequency of fransition is discussed in the next section. 

Nov^' v^t consider the relative fransition probabilities for A^B and A^C. We can define the 
mean frequencies for these specific fransitions as 
TV 

v^^j, = lim and v^^^ = ^ » (15) 

where N^^g and N^^^ the total number of fransitions from state A to state B and state C in the 
time interval [0, r] respectively. Note that since A'^^^ = N^_^g + N^^^., v^_^ = v^_^g + v^_^^ . Then, 
the relative fransition probabilities are given by 



P 



A->C 



^ , (16) 



where P^_^b P a-^c refers to the relative probabilities to the transition A ^ B and A ^ C such that 

Pa-^b+Pa^c =1 • (17) 
Therefore, all the information needed to understand a dynamical system undergoing infrequent events can 
be obtained from the mean fransition frequencies. 

B. Transition state theory and dynamical corrections 

In principle, to obtain the mean frequencies of transition in Eq. (9) and (15), we need to observe a 
statistical number of transitions from the state in question so that we can count the total number of 
transitions to each neighbor and measure the waiting time. Practically it is hard to observe even a single 
transition if we use a conventional dynamics scheme. When two states, usually referred as the reactant 
state and product state, are known, transition state theory can be used to calculate the mean fransition 
frequency between these two states. In this section we briefly review TST fransition rates and the 
dynamical corrections to these rates. 

Given two metastable states A and B, we calculate the mean frequency v^_^g of the fransition 

from A to B. First, we construct a dividing surface S between the corresponding state volimies, which is 
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neither the boundary of A (2/4) nor B (dB) (see Fig. 1), then the mean frequency Vg of crossing S from 
A to B is given by ^ 



where v„ is the velocity normal to the dividing surface S and the prefactor 1/2 is infroduced to account 
for the transition A^B only. Note that Vg > v^^g because the system can recross the dividing surface 
several times before making a transition. No transition has occurred if the trajectory leaving dA re-enters 
A vwthout visiting B. Although v^^e can find the optimized dividing surface vv^hich minimizes Vg 
(variational TST '^), there is still a chance of such an event. 

To account for the possible recrossing of S the dynamical correction factor fj has been 
introduced, ^' ^' which is defined as 



Although we cannot directly calculate u^^g , we can calculate as follows. ^ First, we sample a large 

number of points on the dividing surface S and initiahze velocities according to the correct phase space 
distribution. We initiate a pair of frajectories starting from the same point on the dividing surface, one of 
which has the chosen outward velocity and the other of which has the opposite inward velocity. Then, the 
ratio of the number of trajectories whose outward pair enters B and whose inward pair enters A without 
returning to the dividing surface to the total number of trajectories gives the approximate dynamical 
correction factor. 

C. Hyperdynamics 

Although the original hyperdynamics formulation is based on the TST transition rate, the actual 
fransition rate or the mean frequency of transition does not depend on the choice of the dividing surface. 
In the previous section a dividing surface is used to calculate Og and , but this is merely a construction, 

and the product of these must not depend on the choice of a dividing surface. As we will show it is 
possible to formulate the hyperdynamics method without using a TST dividing surface. In the 
hyperdynamics method it is not necessary to directly calculate the fransition rate. Instead, we need only to 
be able to accurately compute the ratio of the fransition rates in two different potentials and this can be 
done without resorting to the construction of a dividing surface. 




(18) 



u 



(19) 
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where 



First, by rearranging Eq. (9) we have 

N N N 

o,^=hm^ = Mi^^-^ = o,,^xf,^ , (20) 

J— >oo 7- r->oo 7- /V 



v,^^=^^^ , (21) 



/,^=lim^^ , (22) 



dA- 



where iVg^^ is the total number of trajectories crossing dA in the time interval [0, r]. Dg^^ is expressed 



as 

_1 Je^^^J""'""!'^ 



1 L^^J^^I^"!^ 

"-- = 2 Z • ^''^ 



Using N^_^g and iV^^c » /^-> expressed as 

~ Ia-^b Ia-^c (24) 

where 

/,^, = lim ^ and f,_^, = lim . (25) 



Then, we have 



A-^C 



(26) 

'A^B y^"' A^B J J A-^B 

Note that /^^^ and /^^^ depend only on the potential energy values at fi^and in the transition 
region, which is outside of A. To prove this, imagine the distribution of crossing points of the trajectories 
from inward with dA . Since these points also belong to the phase space, they must satisfy the following 
distribution 

Ps(r) = — , (27) 

where = e'^^ dS . Moreover, whether a specific trajectory will return to dA or enter other states 

after it leaves dA is completely determined by the potential energy values in fransition regions. Thus, if 
two potentials have the same values in these regions ( dA and the fransition region), /^^^ and f^^c 
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remain unchanged and the relative transition probabihties in these two potentials are identical. Note that 
although with these two potentials iVg^^ changes, /^^^ , the ratio of iV^.^.^ to iVg^^ , remains 

unchanged. This is the fundamental basis of the hyperdynamics method. 

Now we assume that for a given potential V, we can construct a biased potential Vb such that 
ViP) = V,(f) + AV , (28) 

where a bias potential zlF satisfies the following condition 
f > , in A 

AV(f) = \ . (29) 

[ = , along dA and outside of A 

Then, the ratio of the transition rates in Fand Fj, the boost factor a, is given by 
a = ^ = X= fromEq.(8) 

Ra {tA)b iPA)bl(PA-^)b 



Pa/^sa-^/a-^ 



(PA)bK^dA^)b(fA^)b 

PaI'Jsa^ 

(PA)bKOgA^)b 

f ^dfe-^" 



I ^dre' 

J A+A 



fi-om Eq. (20) 

fi-om /^^ =(fA^)b 

fi-om Eq. (13), (14) and (23). (30) 



Hereafter the subscript b means that the property is obtained in the biased potential. 

It is apparent that the average waiting time in the original potential can be recovered from the 
average waiting time in the biased potential as shown in Eq. (30). Now we consider the probability 
density function of the waiting time. We define the recovered time 7] as ax(tj)i^ . The probability 
density fimction of ;7is given by 

dt] 

a drj a 

= - X (i? J, exp {{R^ ), X (r J J from Eq. (5) for (^ J, 

a 

= e?g)( R^ tj) from Eq. (30) . (31) 
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Therefore, the stochastic outcome of the waiting time in the original potential can exactly be replaced by 
the recovered time tj(= ax(t^)^). 

The state volume A can be well-defined, but A may not. Thus, we derive an alternative formula. 
First, we define the biased boost factor as 

= == • (32) 



Then, using the same derivation in Eq. (30) we have 
f d 

J A 
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a, = . (33) 



r e 



We can also show that t^g= (tAo^t ■ Then, the boost factor becomes 



tA tA,i+tA,o (^bX(tA,i)b+(tA,o)t 

a = - - 



(tA)b (t^,i)b+(tA,o)b (tAj)b+(tA,o)b 
= iPA,i)b^(^b+(PA,o)b > (34) 

where 

(FT \dre-^'^ 

(PA,i)b=== - 7-^ =7 ,^ ' 



itA,i)b+itA,o)b \,~d?e 

J A+A 

Tt ^ \^dre' 

\'-A,o)b J A 



(PA,o)b = 7 — ^"'"r — ^ " — 7^ — W • (^^^^ 
(tA,i)b+(tA,o)b i ^dre-^"^^ 

J A+A 

In the actual hyperdynamics simulations (/)^,.)jand iPA,o)b^^ be approximated by {p^.)l and 
(Pa oYb ^ calculated from the time of residence inside and outside state A in the boosted simulation, i.e. 



Then, we can show that t^=ax (t^ )/, = x (t^ . )/^ + (t^ ^ )^ . Thus, without knowing the actual boost 
factor a we can recover the original waiting time with the biased boost factor . 
In the same way in Eq. (34) we can also show that 
1 1 

a = < . (38) 

PA,i'<^b+ Pa,o Pa,o 
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In the case that p^j , p^ ^ , and a are known, (p^ .)i^md {p^^)^cm be predicted using the following 
relations. 

\ ^dr e \ dr e \ ^dr e 

J A+A J A J A+A 



J A+A J A J A 



Ure-'"! ^dfe-^'^"^"-' ' 

Note that although p^ ^ » p^ ^ , with aggressive boosting (p^ o)iCan be comparable to (p^ . 

Finally, the ensemble averages of an observable 0(r) in the original potential and the biased 
potential have the foUovwng relation 

(0(?)) = -(0(F)e^^^^)^ , (41) 



a 

where 



f ^f(r)e-^''dP 

7 ^ (42) 

^ e dr 

3 A+A 

\ J\r)e''''''dr 



Therefore, the ensemble average of an observable in the original potential can be obtained from the 
ensemble average in the biased potential using Eq. (41) without performing a simulation with the original 
potential. 

III. CRITICAL ISSUES OF HYPERDYNAMICS 
A. Boost factor 



According to Voter's original prescription ^, the time t in the original potential can be recovered 
from the time elapsed in the potential modified with a bias potential A Fusing 

Ntot 

t- Z^^tMD^ , (44) 

i=\ 
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where Ntot is the total number of MD steps within one transition event, is the time interval for 

numerical integration, and is the time at /th MD step ( = /x A?^^). We can prove Eq. (44) using the 
recovered time ( = x a ). By rearranging Eq. (30) we obtain the following relation 

f ^ e-^^ dr f ^ e^^^'^e-^^* dr 

JA+A _ JA+A 

f . e-^"' d? f . e-^"' d7 

JA+A JA+A 



cc= T~' , =^ — = (^^^'"). . (45) 



JA+A JA+A 

From Eq. (45), the boost factor can be approximated by 

e 



^ ^+/?AK[?(?,.)] 



a= e-z'-M . (46) 

TOT 

Note that the MD simulation should be performed with the biased potential, i.e., the force vector should 
be obtained from the biased potential, in order to calculate a using Eq. (46). Since the recovered time can 

be regarded as the original waiting time, we have 
? = 7 = X a 

= Nj^oT X ^t^D X « from = N^oT X ^t^D 

NjOT 

= i:^t^ne'' from Eq. (46). (47) 

(=1 

This agrees with the supposition in Eq. (44). 

In the above derivation it is assumed that Eq. (46) holds, an assumption that largely depends on 
the total number of MD steps sampled in each escape event. If a bias potential is chosen too aggressively, 
the system will stay in the starting state for too short a time to obtain the boost factor a accurately. With 
an inaccurate boost factor, the recovered time loses its statistical meaning and its ability to approximate 
the original time. In an optimistic view expressed in Ref 9, even with the aggressive choice, the 
accumulated time error after many fransitions may vanish because the time error in each transition is not 
correlated with others. However, it is also likely that a bad choice of a bias potential can cause the time 
error in a biased way such that it always yields shorter or longer estimates at every fransition. Moreover, 
using too conservative a bias potential is not desirable for obvious reasons. 

Since we need the accurate boost factor for each transition, a better sampling scheme for the 
boost factor is needed and this sampling may be performed as a pre-simulation before the actual 
simulation. By rearranging Eq. (30), we have 
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f . e-P'' dr \ ^ e-^^ dr \ . e"^^ dr 

^ _ JA+A _ JA+A JA+A 

f ^ e-^"' dr f . e-^^' ti? f . e"^^ c/F 
J^+^ J^+^ J^+^ 

f e-^'^ dr f e^^^^-^'V^^ ^/r 

_ JA+A JA+A 

l-fi(_W-V) \ 



^ , (48) 



where W is an intermediate bias potential. Thus, in principle we can use any type of intermediate bias 
potential (even the original potential) to calculate the boost factor. 

Let us consider a case with a bias potential using a local variable A (e.g. the lowest eigenvalue of 
the Hessian matrix). From Eq. (48) we have the following relation 

a = i = i (49) 

J —GO 

where p (A,) is the probability density function of /I in A + A with the original potential. If the state 
volume is defined by A < A^^ (i.e. AV = when A > A^^), then we have 

a = ^- ^- < ^- . (50) 

r /7(/l)e-^^^Wj^ f -p(A)e-^^'^'^dA + piA>A) Pi^> Kr) 

J —00 J —GO 

Therefore, the boost factor cannot exceed \l p(^X> A^^~) . If we use an intermediate bias potential 
W = V + AW (A) using the same local variable A for the sampling, the probability density flmction 
(A) in Wis given by 

l^./iAir)-A,)e-^-df £^/(A(r-) - ^)e-^(-™JF 



Pw (■^o) = 



JA+A _ JA+A 

f .e-'^^'dr f ^e-'^^'dr 

JA+A JA+A 

e-^^""''"' x^^^^SiAif)-Ao)e-^''dr 
f .e-^^'dr 

JA+A 

g^wio)xf ^SiAiP)-Ao)e-^''dP \ ^e-^'dP 
f .e-^Vr f ^e-^^'dr 

JA+A JA+A 
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where 



<^w=^f^^. • (52) 



iA+A 

3hW 



Due to the prefactor e in Eq. (51), p^i?^ is biased such that the region in the configuration space 

where AfFis larger is less sampled, and the region having smaller AfFis more sampled. Thus, we can 
effectively sample the configuration space choosing a proper intermediate bias potential W. For some 
cases, more than two bias potentials can be used for sampling and the resultant probability density 
functions can be combined, for example, using the weighted histogram analysis method (WHAM) as in 
umbrella sampling. Note that once we construct the probability density function in one specific potential, 
the boost factors in any bias potentials can be calculated. When using pre-simulations to calculate the 
boost factor, it is critical to have sufficient numbers of MD steps for sampling in order to obtain an 
accurate boost factor. 

B. State volumes and bias potentials 

Since the original formulation of hyperdynamics '^used the TST transition rate, the bias potential 
was constructed based on the TST dividing surface. The conventional TST dividing surface uses the 
steepest ascent/descent path described by 

^^ = ±^ , (53) 

ds |VF| 

where ds (=\dr \) is the arc length of the iA'-dimensional curve in the configuration space. All the 
points which can be led to a minimum by the steepest descent path comprise the state defined by the 
minimum, and the points on the boundary, which converge to one of the first-order saddle points instead 
of the minima, define the TST dividing surface. Instead of this conventional TST dividing surface Voter 
used an approximate dividing surface defined by 

Ci-VK = and < , (54) 

where is the lowest eigenvalue of the Hessian matrix and Cj is the corresponding eigenvector. 

However, calculating the eigenvalue, the eigenvector, and their derivatives (to obtain force) is extremely 

time-consuming. As stated in the previous sections, our hyperdynamics formulation does not use a TST 

dividing surface, but depends on the boundaries of the state volumes. This gives more flexibility in how 
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we define bias potentials. To illusti-ate this point we will next construct a condition for a state volume and 
a corresponding bias potential. 

The requirements for defining a state volume are that it contains a region where the probability 
densities are concentrated and that it must be distinguishable from other such volumes, insofar that it must 
not share any position vector with these other states. It is desirable that the system exits the state volume 
only when it makes a tiansition in order to maximize the boost factor as shown in Eq. (50), but this 
property is not required and does not affect the accuracy of the hyperdynamics simulation. 

If a local variable satisfying the above requirements can be found, then we can construct a bias 
potential as a function of this local variable. For example, as in Voter's original bias potential, we can use 
the lowest eigenvalue s^of the Hessian because each region satisfying > Oin the configuration space 

satisfies these requirements. However, even the most efficient proposed schemes for the eigenvalue 
calculation require multiple computations of a force vector, and are therefore very computationally 
intensive. '"'^^ 

One possible alternative choice is the potential energy slope or curvature along the direction 
vector connecting a configuration r and the potential energy minimum . They are defined by 

dV 

s^VV-u= , (55) 

d u 

c = —-r , (56) 



du' 



where 




(57) 



3N 



Y^^h-ro^kf . (58) 

where is a component of the position vector r and ^ is a component of the position vector at 

the minimum. In some specific systems, the slope in Eq. (55) and/or the curvature in Eq. (56) become 
smaller when the system approaches a low probability density region (near the boundary of a state volume) 
and in such cases it may be possible to determine critical values for them to define the boundary of the 
state volume. The derivatives are given by 

d^ V gi'U 



87 I 



du^ I 



(59) 
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dr I 



du^ I 



(60) 



where 



du du 



d\Fl) d'F dF 

g2= V^ = --r^^-2— > (62) 

du du du 



and 



F = -VV . (63) 
Note that all the higher order derivatives along a specific direction can easily be approximated using a 
finite difference scheme. 

A more inexpensive choice is the hyper-distance irom the minimum defined in Eq. (58). The 
derivative is given by 

dl r -r„ 
or I 

Note that the extra computations for this derivative are negligible. If we observe the hyper-distances 
and Ig from two minima A and B respectively and the probability distributions of I ^ and Ig are 
distinguishable depending on whether the system stays at A or B, then this hyper-distance can be used to 
construct a bias potential. However, we should approach this method with caution because in some cases 
the probability distribution of the hyper-distance in one state may not be distinguishable from that in 
another state, hi general the thermal fluctuation of the hyper-distance from a minimum will rise as the 
number of atoms in the system increases and in some systems transitions occur in a localized region 
including a relatively small number of atoms. In such cases the hyperdistance between two minima 
depends only on these few atoms and the hyperdistance associated vwth the transition is likely to be 
smaller than the magnitude of the combined thermal fluctuations in the system. However, even in this 
case the hyper-distance can be used if we can identify where a transition will occur. For example, in the 
firictional sliding system with a sharp tip scanning a substrate, important transitions always occur at the 
interface and we can use the atoms at the interface to measure the hyper-distance rather than including all 
atoms. 



IV. APPLICATION 
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A. Model 

The methodology described above has been tested with a 3-dimenstional model of an Atomic 
Force Microscope with a tip and a substrate as shown in Fig. 2 (a). The substrate is a FCC crystal 
consisting of 1800 atoms and the tip consisting of 183 atoms is made by carving a FCC crystal into a 
conical shape with flat ends. The tip has the same lattice parameter as the substrate and is joined to the 
substrate in the [001] direction. As shown in Fig. 2 (b) nine atoms on the bottom layer of the tip are in 
contact with the substrate. Since the tip and the substrate have the same lattice parameter and are aligned 
along the same orientation, the tip atoms are in registry with the substrate. All the quantities are expressed 
in the length units of a; the energy units of s, and the mass units of m. Time is measured in units of 
r = o" (w / f')"^ . Hereafter the units are omitted unless an ambiguity would arise. 

The atoms on the bottom layer of the substrate are fixed to prevent rigid body translation, but they 
interact with other atoms. Periodic boundary conditions are apphed in the horizontal directions. The 
relative motions of the atoms on the top layer of the tip are constrained, but they can move like a rigid 
body. The stiffness of an AFM cantilever is modeled by linking a spring with a stiffness of 10 to the 
top layer of the tip and the spring is elongated by changing the position of a shder located at the end of 
the spring. A normal force of 5 is also applied downward as shown in Fig. 2 (a). 

The interactions of the atoms are modeled by the Lennard- Jones potential, 

F(r) = 4f, 

where Sab is the bond energy between the atom of the type a and the atom of the type h, Uah is the 
characteristic length parameter, and r is the distance between the two atoms. We used the following 
length scale and energy parameters. 

^ss ^^tt ^^ts = 1 -0 ' = =1.0, s^^ =0.2 (s: substrate, t: tip) 
Note that we used a smaller bond energy for the interaction between the tip and the substrate to guarantee 
that the rearrangement of atom positions always occurs at the interface rather than inside the tip. We 
observe the dynamics of this system while it is maintained at a constant temperature (7'= 0.01 ^ks) by a 
Nose-Hoover chain thermostat. The equations of motion are solved using a modified velocity-Verlet 
algorithm. 

In a sliding simulation modeling an AFM experiment the slider moves at a constant velocity, 

but in this study we used a fixed slider to test the methodology with a non-driven system with constant 

boundary conditions. The extension of the current methodology to a driven system can be found in Ref 
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24. As shown in Fig. 2 (b) the tip can hop to one of eight neighboring positions, each of which defines a 
metastable state. We biased the relative transition probabilities to these neighboring states by elongating 
the spring in the positive x direction (xs = 1.3) so that the hop in the positive x direction dominates other 
hops. Hereafter the initial state is referred to state A and the final state in the positive x direction is 
referred to B. We observed a transition fi"om A to B by performing MD simulations with the original 
potential (conventional simulations) and the biased potentials (hyperdynamics simulations). A transition 
fi-om state A to state B is detected by performing periodic minimizations using a scheme called FIRE. 

From the pre-simulation we have observed that whether the system is at state A or state B can be 
distinguished by measuring the hyper-distances fi-om each minimum. Therefore, we constructed the bias 
potentials using the hyper-distance and the following functional form 

2 

l^<l <l^j , (66) 
l>l^ 

where = 0.2 , = 1.0 , and = 4. We prepared four different bias potentials (r/°' , V^"" , r/"" , 
r/°^), each of which has AF^^^ = 0.02, 0.04, 0.06, 0.08 respectively. For the original potential and 
four biased potentials, we simulated 100 samples with different initial conditions. 

B. Results 



^v{l) = 



1- 



First, we discuss the results of boost factor calculation. The metastable state volume A is defined 
as a hyper-sphere in the configuration space, centered at the initial minimum and with a radius of 1^ . If 

the probability density fiinction yO^ (/) of the hyper-distance / in a biased potential W(= V + Afr(/))is 

known, then the biased boost factor of any other biased potential V^^{=V + AF(/)) , defined in Eq. (32) 
and Eq. (33), is given by 

«6 = —, , (67) 

J 

where (I) is normalized such that 
\l"pAl)dl = l . 
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For a thermodynamic sampling to obtain (/) in a biased potential Wwe performed a MD simulation 

with W. Moreover, by definition, we must not observe any transition during the sampling. Thus, 
configurations are periodically stored at prescribed time steps and when a transition is detected we resume 

the simulation staring fi'om the most recently stored configuration after assigning a new set of Boltzmann- 
distributed velocities. 

Fig. 3 shows the probability distributions of the hyper-distance from the minimum A in the 
original potential obtained by the MD simulations consisting of 1,000,000 and 5,000,000 integration steps 
respectively. As the number of MD steps increases, the distribution becomes more accurate. However, 
even with the larger number of steps (5,000,000), the distribution in the transition region A (I > l^) is 
noisy because this region is rarely visited during the sampling with the original potential. 

In Fig. 4, the distributions in various biased potentials (V°°^ ,V°°\V°°\V°°^) are compared 
with the distribution in the original potential. All the data were obtained by the MD simulations with 
5,000,000 steps except for the distribution in F/°** , which was obtained after 3,000,000 steps. With an 

aggressive biased potential like °* , we observed many transitions during the sampling. At these times 
the sampling had to be resumed from the previously stored configuration resulting in a significantly 
increased simulation time. One thing to note in these graphs is that as A F^^,^ increases, the probability to 

be in A increases. The role of a biased potential is to force the system to move into this region, increasing 
the probability for a transition. Moreover, eventually with V^ °^ we have two peaks in the distribution. It 

is desirable to avoid the creation of such an extra peak because if these two peak regions are separated by 
very low probability regions, each peak becomes a separate state which can impede ergodically sampling 
the phase space. 

Once we have the probability distribution in one of these potentials (V , V°-°^, V°-'^, V°°^ , V°-'^), 
we can calculate the probability distribution and the boost factor in any other potential. Fig. 5 shows the 
biased boost factors in F/ , F^*""* , V^ °^ , F/ as functions of integration steps, and each graph has two 
curves: one is obtained from the original potential and the other fi-om the corresponding biased potential. 
For the same reason above, the biased boost factor in F/ °* calculated from the distribution in V°"^ was 

obtained only up to 3,000,000 steps. In principle these two boost factors from the original potential and 
the biased potential should converge as the number of steps increases, but as shown in Fig. 5 the 

difference becomes large as A F^^^^ increases. This can be explained by referring to Fig. 6, which shows 

the probability distributions of the hyper-distance in V^'°^ recovered from the original potential and the 

19 



biased potentials. If we obtained the exact distribution in each sampling with a specific potential, the 
recovered distributions in other potentials should also be exact. However, because we use a finite number 
of steps for a sampling, there is always an inaccuracy in the distribution. This inaccuracy also makes the 
error in the boost factor calculation larger as A increases because of the exponential term in boost 

factor calculation. 

As an alternative sampling, we used the umbrella sampling method. In the umbrella sampling 
instead of using either the original potential or one of the biased potentials for the hyperdynamics 
simulations, we used harmonic potentials given by 



where k' =k^ =k^ =k^ =k' =1.0 , and =0.4 , =0.6 , =0.8 , =1.0 , =1.2 . Each 
harmonic potential is designed to sample the region near its center more accurately. For each harmonic 
potential we sampled 1,000,000 points and once we have the probability distribution of the hyper-distance 
in each harmonic potential, we can reconstruct the distribution in the original potential using WHAM as 
shown in Fig. 7. Although we used the same number of MD steps both in the conventional sampling 
and in the umbrella sampling (5x 1,000,000 = 5,000,000), we obtained a more accurate result with the 
umbrella sampling. 

Fig. 8 shows the biased boost factors of the hyperdynamics biased potentials ( F/"^ , V^'''* , V^'''^ , 
yom ^ jjjg boost factors calculated from the conventional sampling with the original potential and the 
hyperdynamics biased potentials are compared with the boost factor obtained from the umbrella sampling. 
The results with significantly deviate from the imibrella sampling boost factors because of the 

reasons explained above. 

Next, we performed the MD simulations to see the transition from A to B with the original 
potential and the biased potentials. For each potential we prepared 100 samples with different initial 
conditions and measured the waiting times. Fig. 9 (a) and (b) show the average waiting time in the state 

volume (/</[/), i , and the average waiting time in the transition region (l > l^), t^ ^. As AF^^,^ 

increases, that is, the boost factor increases, ^ decreases, but ^ remains almost within the error bar, 
which is the standard error, regardless of the potentials as we discussed in Sec. II. C. Fig. 10 shows the 



(/ = 1,2,3,4,5) 



(68) 
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average number of crossings of dA(l = l^j), the boundary of the state volume, before a transition. It also 

remains within the standard error, which confirms f^^^ — (fA^B)h ■ 

Finally, Fig. 1 1 shows the original waiting times recovered from the waiting times measured in 
the hyperdynamics simulations with the biased potentials through two different methods. In one method 
we use the boost factor obtained from the pre-simulation as discussed in Sec. III. A, the original waiting 
time is recovered using the relation, cc^ x (t^ + (t^ „)j , and the other method by Voter uses Eq. (44), 

that is, the boost factor is calculated during the simulation on the fly. The original waiting time is well 
recovered from the waiting times in V^'"^ , V°''^* , V^'°^ , but the recovered waiting time from the 

simulations with F/°* deviates with a larger error even with the pre-calculated boost factor. We 

conjecture that it is because the biased waiting time in V^'"^ becomes too short to satisfy the criterion that 

the system is fully equilibrated in a state before making a fransition. The waiting times recovered using 
Eq. (44) show larger deviations from the results obtained from the pre-calculated boost factor. 

V. CONCLUSION 

We have reformulated the hyperdynamics method, which was devised to accelerate the 
conventional MD simulation, in a rigorous way that does not require a TST dividing surface. First, we 
define a transition as an event between two state volumes in the configuration space rather than a crossing 
of a TST dividing surface. By the fact that the ratio of the number of transitions to the number of 
crossings of the boundary of these state volumes remains unchanged in a biased potential, we showed that 
the boost factor, which is the ratio of the transition rates in the biased potential and in the original 
potential, can be used to exactly recover the original waiting time. 

We presented a new perspective to see the boost factor as a multiphcation factor between the 
waiting times in the biased potential and in the original potential so that it can be calculated in a pre- 
simulation rather than during the simulation as in Voter's original method. To accurately calculate the 
boost factor we discussed the various aspects of thermodynamic sampling and concluded that the 
imibrella sampling gives a more accurate result when compared with the conventional sampling 
performed "on-the-fly". 

Moreover, with the criteria of a state volume rather than using a TST dividing surface, we devised 

new bias potentials using local variables other than the eigenvalue and the eigenvector of the Hessian 

mafrix. Among these local variables the hyper-distance from a minimum turned out to be the most 

efficient. However, an important caveat must be noted. The hyper-distance can only be used in a system 
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where the hyper-distances between minima are distinguishable from the thermal fluctuation of the hyper- 
distance from each minimum, a condition that holds for our AFM simulations because a confined 
fransition region can be identified a priori. The new methodology was tested for an AFM system 
modeled by the Lennard- Jones potential, and the simulation gave the results consistent with the theory. 
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FIG. 1. An example of 2-dimensional potential energy contours (thick solid curves). The 
entangled light solid curve represents a system trajectory, the dashed circles represent the 
boundaries of state volumes, and the dotted curve is a TST dividing surface. Note that a 
transition from one metastable state to another is evident although several crossings of the 
boundaries of these state volumes occur before the system makes a transition. 
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(a) (b) 

FIG. 2. A diagram of 3 -dimensional AFM model consisting of a tip and a substrate, (a) The 
normal force and the spring force are applied on the top layer of the tip and the spring force 
can be controlled by changing the slider position x.,. (b) The configuration of atoms at the 
interface between the bottom layer of the tip and the top layer of the substrate. The arrow 
shows the direction for a transition to one of the neighboring states. 
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FIG. 3. The probability distribution of the hyper-distance from the minimum A in the 
original potential. The dashed curve was obtained after 1 ,000,000 MD steps and the solid 
curve was obtained after 5,000,000 MD steps. The dashed vertical line shows the boundary 
of the state volume A. 
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FIG. 4. The probability distributions of the hyper-distance from the minimum A in the 
biased potentials ( F^" , V°°'* , V°°^ , F," ). For comparison, the distribution in the original 

potential is plotted together with the distribution in the biased potentials. All the results 
were obtained by the MD simulations with 5,000,000 steps except for the distribution in 



F/"' , which was obtained after 3,000,000 steps. 
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FIG. 5. The biased boost factors in the biased potentials ( °% , T," "' ). For 

comparison, the biased boost factor calculated from the hyper-distance distribution in the 
original potential is plotted together with the biased boost factor obtained from the hyper- 
distance distributions in the corresponding biased potentials. 
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FIG. 6. The probability distributions of the hyper-distance in F," "** . One was obtained from 
the simulation with F^"** and the others were recovered from the distributions in the 
original potential and F," , F,""" , F," . 
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FIG. 7. The probability distribution of the hyper-distance in the original potential obtained 
from the conventional sampling and the umbrella sampling using 5 harmonic potentials. 
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FIG. 8. The biased boost factors of the biased potentials ( V'-°' , V°-°' , V° °' , V° '' ). The 
open dots show the boost factors obtained from the conventional sampling with the 
potential whose AF^^^^is given in the x axis and the sohd dashed lines show the boost 
factors obtained from the umbrella sampling. 
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FIG. 9. The average waiting times in the original potential and the biased potential whose 
AFj^^j^is given in the x axis, (a) The average waiting time in the state volume A, . . (b) 

The average waiting time in the transition region A , ^ . The red dashed line is the total 
average and the error bars show the standard error. 
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FIG. 10. The average number of crossing of the boundary of the state volume A (dA) in 
the original potential and the biased potential whose AF,^^^is given in the x axis. The red 
dashed line is the total average and the error bars show the standard error. 
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FIG. 1 1 . The average of the original waiting times recovered from the waiting times in the 
biased potentials, whose AF^^^is given in the x axis, with the boost factors calculated via 

two different schemes; one is the pre-simulation and the other is the on-the-fly calculation. 
The red dashed line is the average waiting time in the original potential and the error bars 
show the standard error. 
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